Gulf of Maine Current Year SST Report

Temperature tracking for the Gulf of Maine

Author
Affiliation

Gulf of Maine Research Institute

Published

August 29, 2022

Gulf of Maine 2022: Sea Surface Temperature

This report was created to track the sea surface temperature regimes for marine regions of interest to the Gulf of Maine Research Institute. Historically the focus has been on a central snapshot of the Gulf of Maine. A region that has experienced profound temperature increases in recent years.

Satellite sea surface temperature data used was obtained from the National Center for Environmental Information (NCEI). With all maps and figures displaying NOAA’s Optimum Interpolation Sea Surface Temperature Data.

DISCLAIMER: Any data within 2-weeks of the current date are subject to revision and may change. Please use caution when reporting information that contains these values.

Gulf of Maine On a Map

Whenever discussing the Gulf of Maine in this report, we refer to the following spatial extent displayed below. The coordinates of this bounding box are the same coordinates used to clip the sea surface temperature data.

Code 
# Load the bounding box for Andy's GOM to show they align
poly_path     <- region_paths[[params$region]][["shape_path"]]
region_extent <- st_read(poly_path, quiet = TRUE)

# Pull extents for the region to set crop extent
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]


# Map the study area
map_study_area(region_extent = region_extent,
               x_buffer =  c(-2, 2),
               y_buffer =  c(-0.75, 0.75), 
               new_england_sf = new_england, 
               canada_sf = canada,
               greenland_sf = greenland)

1: Temperature History

Using the coordinates shown above we can create an area-specific temperature history. For any day of the year since September of 1981, data is available to calculate the average sea surface temperature within that area for each day.

This report compares observed temperatures against the expected conditions based on a 30-year climatology. The standard reference period used for the climatology here is 1982-2011, as used here.

Code 
# Load the timeseries
timeseries_path <- region_paths[[params$region]][["timeseries_path"]]
region_timeseries <- read_csv(timeseries_path, 
                              col_types = cols(), 
                              guess_max = 1e6) %>% 
  mutate(tod = format(time, format = "%H:%M:%S")) 


# Clean up the data - add labels
region_timeseries <- region_timeseries %>% 
  mutate(time = as.Date(time)) %>% 
  distinct(time, .keep_all = T) %>% 
  supplement_season_info() %>% 
  filter(year %in% c(1982:2022))



# Summarize by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>% 
  filter(year %in% c(1982:2021)) %>% 
  group_by(year) %>% 
  temperature_summaries() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")),
         anom_direction = ifelse(area_wtd_anom > 0, "Hot", "Cold"))


# Get heatwave statuses for each day:
# Uses area weighted sst by default
# Set up the date within a year for Heatmap

base_date <- as.Date("2000-01-01")
region_hw <- pull_heatwave_events(
  temperature_timeseries = region_timeseries, 
  threshold = 90, 
  clim_ref_period = c("1982-01-01", "2011-12-31")) %>% 
  supplement_hw_data() %>% 
  filter(doy != 366) %>% 
  mutate(year = year(time),
         yday = yday(time),
         yr_rev = factor(year),
         yr_rev = fct_rev(yr_rev),
         flat_date = as.Date("2000-01-01") + yday - 1)



 # Current Calendar Year
this_yr_hw <- region_hw %>% 
  filter(year(time) == 2022)

# Hottest year data
hottest_yr <- filter(region_hw, year(time) == 2021)


# Global Temperature Anomaly Rates
global_anoms <- read_csv(
    paste0(oisst_path, "global_timeseries/global_anoms_1982to2011.csv"), 
    guess_max = 1e6,
    col_types = cols()) %>% 
  mutate(year = year(time)) %>% 
  filter(between(year, 1982, 2021))

# summarize by year again
global_summary <- global_anoms %>% 
  group_by(year) %>% 
  temperature_summaries() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")),
         anom_direction = ifelse(area_wtd_anom > 0, "Hot", "Cold"))

Climatological Patterns of the Gulf of Maine

Daily climatologies are a record of daily averages based on the day of the year. They record what the average temperature has been for each day of the year, across the specified reference period.

In addition to the daily average, we also look at how variable temperatures are. This variability is used to benchmark how “rare” extreme events are to determine whether they are part of a natural cycle or not.

Common benchmarks for extreme events are the 10th and 90th percentile. Temperatures above the 90th percentile are among the hottest 10% of days in the reference period. Temperatures below the 10th percentile correspond to the coldest 10%. The remaining 80% of days fall somewhere in-between and showcase the range of temperatures we might expect to occur given the natural variability of the climate.

Code 
# Set plot colors:
clim_colors <- c(
  "Daily Temperature"      = "gray80",
  "Climatological Average" = "gray10",
  "90th Percentile"        = "coral3",
  "10th Percentile"        = "royalblue")

# pull one year for the lines, the reference period for points
one_yr <- region_hw %>% filter(yr == 2021)
ref_yrs <- region_hw %>% filter(yr %in% c(1982:2011))

# Make a plot
ggplot(data = one_yr, aes(x = flat_date)) +
    geom_line(data = ref_yrs, aes(y = sst, group = year, color = "Daily Temperature"), alpha = 0.6) +
    geom_line(aes(y = seas, color = "Climatological Average"), size = 1) +
    geom_textline(aes(y = mhw_thresh, color = "90th Percentile", label = "Heatwave Threshold"), linetype = 1, size = 4, hjust = 0.4, linewidth = 0.5) +
    geom_textline(aes(y = mcs_thresh, color = "10th Percentile", label = "Coldspell Threshold"), linetype = 1, size = 4, hjust = 0.6, linewidth = 0.5) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = expansion(add = c(0,0))) +
  scale_y_continuous(labels = number_format(suffix = " \u00b0C")) +
    scale_color_manual(values = clim_colors) +
    labs(x = "Day of the Year", y = "Sea Surface Temperature",
         title = str_c(tidy_name, " Daily Climatology"),
         subtitle = "Anticipated Daily Temperatures Based on a 1982-2011 Climate") +
    theme_gmri() +
    theme(
      legend.position = "bottom",
      legend.key.height = unit(.5, "lines"),
      legend.key.width = unit(5, "lines"),
      plot.margin = unit(c(.5,1,.3,.5), "cm"),
      plot.title.position = "plot",
      legend.margin = margin(c(7,0,0,0)),
      legend.justification = "center") +
      guides(color = guide_legend(
        title = "",
        title.hjust = 0.5,
        nrow = 1,
        title.position = "left",
        label.position = "top", 
        override.aes = list(
          color = c("gray50", "gray10", "coral3", "royalblue"),
          shape = c(NA, NA, NA, NA),  
          linetype = c(1, 1, 1, 1),
          alpha = c(1,1,1,1), 
          label = c("", "", "", "")),
        byrow = T))

Working in Anomalies

Using the climatology as a reference, we can see hot much “hotter” or “colder” a given day is than what we would on average expect. This difference from the average is what we call anomalies.

In the table below, Sea Surface Temperature is the mean temperature observed for that date averaged across all cells within the area. Climate Avg. & Climate SD are the climate means and standard deviations for a 1982-2011 climatology. Temperature Anomaly is the daily observed sea surface temperature minus the climate mean.

In this way a Temperature Anomaly is just: How much hotter or colder is it than the average temperature for that day.

Code 
# Display Table of last 6 entries
tail(region_timeseries) %>% 
  mutate(across(where(is.numeric), round, 2)) %>% 
  select(
    Date = time,
    `Sea Surface Temperature` = sst,
    #`Area-Weighted SST` = area_wtd_sst,
    `Day of Year` = modified_ordinal_day,
    `Climate Avg.` = sst_clim,
    #`Area-weighted Climate` = area_wtd_clim,
    `Temperature Anomaly` = sst_anom#,
    #`Area-Weighted Anomaly` = area_wtd_anom
    
  ) %>% gt() %>% 
    tab_header(
    title = md(paste0("**", tidy_name, " - Regional Sea Surface Temperature", "**")), 
    subtitle = paste("Temperature Unit: Celsius")) %>%
  tab_source_note(
    source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Climatology Reference Period: 1982-2011.*"))
Gulf of Maine - Regional Sea Surface Temperature
Temperature Unit: Celsius
Date Sea Surface Temperature Day of Year Climate Avg. Temperature Anomaly
2022-08-20 18.38 233 17.03 1.35
2022-08-21 18.98 234 17.00 1.98
2022-08-22 19.16 235 16.98 2.18
2022-08-23 19.48 236 16.89 2.58
2022-08-24 19.63 237 16.76 2.87
2022-08-25 19.87 238 16.70 3.18
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Climatology Reference Period: 1982-2011.
Code 
# march 1st sst
mar1 <- region_timeseries %>% 
  filter(modified_ordinal_day == 61) %>% 
  distinct(sst_clim) %>% 
  pull(sst_clim)

Detecting Long-Term Changes

One reason we go through this trouble of calculating climatologies and anomalies is to set expectations around how much we expect things to naturally vary, and to detect when things fall outside that range.

In the case of the Gulf of Maine we are now consistently outside of the range of temperatures normally expected from 1982-2011. The following plot colors the monthly average by how far temperatures are from the climatological average. Blue colors represent a month that was cooler than average, with red indicating a warmer than average month.

Code 
# Make Monthly Averages
monthly_temps <- region_timeseries %>% 
  group_by(year = year(time),
           month = month(time)) %>% 
  summarise(sst = mean(area_wtd_sst),
            anom = mean(area_wtd_anom),
            .groups = "drop") %>% 
  mutate(time = as.Date(str_c(year, str_pad(month, width = 2, pad = "0", side = "left"), "01", sep = "-")))


# Plot monthly timeline
ggplot(monthly_temps, aes(x = time, y = 0)) +
  geom_tile(aes(fill = anom)) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), oob = oob_squish) +
  scale_y_continuous(expand = expansion(add = c(0,0.01))) +
  scale_x_date(expand = expansion(add = c(45,45))) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
        plot.subtitle = element_text(hjust = 0.125),
        panel.border = element_rect(color = "black", size = 1, fill = "transparent"),
        axis.ticks.y = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        panel.grid = element_blank(),
        plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40")
        ) +
  labs(title = str_c(tidy_name, "\nMonthly Sea Surface Temperature Anomalies"),
       x = "",
       caption = "Monthly Averages from 1982-2021  |  Climatology Period 1982-2011")

Warming Rates

Seeing this change towards a warmer climate, it is then natural to ask how rapidly the change is happening. This is where we turn to warming rates. The warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2021.

The overlaid trend lines then track how warming has increased with time. A dotted line has been included to show how the global average temperature has changed during the same period.

Seasonal Patterns

Code 
# Doing seasons by meteorological Definitions
quarter_summary <- region_timeseries %>% 
 filter(year >= 1982) %>% 
  group_by(year = season_yr, season_eng) %>% 
  summarise(sst = mean(sst, na.rm = T),
            sst_anom = mean(sst_anom, na.rm = T), 
            area_wtd_sst = mean(area_wtd_sst, na.rm = T),
            area_wtd_anom = mean(area_wtd_anom, na.rm = T),
            .groups = "drop") %>% 
  mutate(season_eng = factor(season_eng, levels = c("Winter", "Spring", "Summer", "Fall")))

# Plot
quarter_summary %>% 
  ggplot(aes(year, area_wtd_anom)) +
  geom_line(group = 1, color = "gray60", linetype = 3) +
  geom_point(size = 0.75) +
  geom_smooth(method = "lm", 
              aes(color = "Regional Trend"),
              formula = y ~ x, se = F, linetype = 1) +
  stat_poly_eq(formula = y ~ x,
               color = gmri_cols("orange"),
               aes(label = paste(..eq.label..#, ..rr.label.., sep = "~~~"
                                 )),
               parse = T) +
  scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
  labs(x = "", 
       y = "Sea Surface Temperature Anomaly",
       caption = "Regression coefficients reflect annual change in sea surface temperature.") +
   scale_y_continuous(labels = number_format(suffix = " \u00b0C")) +
  theme_gmri() +
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent")) +
  facet_wrap(~season_eng, ncol = 1)

Overall Temperature Increase

Code 
dat_region <- annual_summary %>% 
  filter(year %in% c(1982, 2021)) 
dat_global <- global_summary %>% 
  filter(year %in% c(1982, 2021)) 
dat_list <- list(dat_region, dat_global) %>% setNames(c(tidy_name, "Global Oceans"))
dat_combined <- bind_rows(dat_list, .id = "Area") %>% 
  select(Area, area_wtd_sst, year) %>% 
  pivot_wider(names_from = year, values_from = area_wtd_sst)


# 
ggplot(dat_combined, aes(x = `1982`, xend = `2021`, y = fct_rev(Area))) +
  geom_dumbbell(colour = "lightblue", 
                colour_xend = gmri_cols("gmri blue"), 
                size = 3, 
                dot_guide = TRUE, 
                dot_guide_size = 0.5) +
  geom_text_repel(
    aes(x = `2021`, y = fct_rev(Area), 
        label = str_c("+", round(`2021`- `1982`, 2), " C")),
                color = gmri_cols("gmri blue"),
                vjust = 4,
                hjust = 0,
                segment.size = 0.5,
                seed = 123) +
  labs(x = "Sea Surface Temperature", 
       title = "Change in Sea Surface Temeprature - 1982-2021", 
       y = "") +
   scale_x_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C"))

Marine Heatwaves

A marine heatwave is defined as a situation when seawater temperatures exceeds a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event. The heatwave threshold used below was 90%. The heatwave history for Gulf of Maine is displayed below:

Note: For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package {heatwaveR}.

Heatwave Events

Code 
# Plot heatmap
heatwave_heatmap <- heatwave_heatmap_plot(hw_dat = region_hw, temp_units = "C", end_yr = 2022)


# Assemble pieces
heatwave_heatmap

2: 2022 Observations

2022 in Context

When 2022 is overlaid against the daily values and climatological average for the region, we can see how temperatures compare to the history of the region.

Code 
# Number of heatwave events this current year
num_hw_days <- sum(this_yr_hw$mhw_event, na.rm = T)


# Plot this year among all previous

# Make a plot
ggplot(data = one_yr, aes(x = flat_date)) +
    geom_point(data = region_hw, aes(y = sst, color = "Daily Temperature"), alpha = 0.6) +
    geom_textline(
      aes(y = seas, 
          label = "Climatological Average", 
          color = "Climatological Average"),
      linetype = 1, size = 4, hjust = 0.5) +
    # Current Year
    geom_line(data = this_yr_hw,  
              aes(x = flat_date, sst, color = "2022"), size = 1) +
    # Hottest Year on Record 
    geom_line(data = hottest_yr,  
              aes(flat_date, sst, color = "2021"), size = 1) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = expansion(add = c(0,0))) +
    scale_color_manual(
      values = c(
        "Daily Temperature" = "gray80",
        "Climatological Average" = "gray10",
        "2021" = "darkred",
        "2022"= as.character(gmri_cols("orange"))
      )
    ) +
    scale_y_continuous(labels = number_format(suffix = " \u00b0C")) +
    labs(x = "Day of the Year", y = "Sea Surface Temperature",
         title = str_c(current_yr, " Among the Hottest for the Gulf of Maine"),
         subtitle = "Conditions suggest a persistant warm regime, comparable to the Gulf of Maine's hottest year, 2021.") +
    #theme_gmri() +
    theme(
      legend.position = "bottom",
      legend.key.height = unit(.5, "lines"),
      legend.key.width = unit(5, "lines"),
      plot.margin = unit(c(.5,1,.3,.5), "cm"),
      plot.title.position = "plot",
      legend.margin = margin(c(7,0,0,0)),
      legend.justification = "center") +
      guides(color = guide_legend(
        title = "",
        title.hjust = 0.5,
        nrow = 1,
        title.position = "left",
        label.position = "top",
        override.aes = list(
          shape = c(16, NA, NA, NA),
          linetype = c(0, 1, 1, 1),
          alpha = c(1, 1, 1, 1),
          label = c("", "", "", "")),
        byrow = T))

Code 
# Seasonal Comparisons:

# Compared to 2022:
year_comp <- region_hw %>% 
  split(.$season_yr) %>% 
  map_dfr(function(x){
  
  # Make sure you only use the available days from this year
  x_filtered <- filter(x, yday(time) <= yday(max(this_yr_hw$time)))
  
  # Check each season
  x_filtered %>% 
    split(.$season_eng) %>% 
    map_dfr(function(.y){
      .y %>% 
        summarise(
          `Average Temperature` = mean(sst),
          `Temperature Anomaly` = mean(sst_anom),
          `Seasonal Low` = min(sst),
          `Seasonal High` = max(sst),
          `Days` = n(),
          `HW Days` = sum(mhw_event)) %>% 
        mutate(`Percent HW` = (`HW Days`/`Days`)*100) %>% 
        mutate(across(where(is.numeric), .fns = round, 2)) %>% 
        select(-c(`Days`, `HW Days`))
          
    }, .id = "Season")
  
  
}, .id = "Year") 


# Do some arranging
year_comp %>% 
  group_by(Season) %>% 
  arrange(desc(`Average Temperature`))  %>% 
  mutate(Rank = row_number()) %>% 
  ungroup() %>% 
  filter(Year == 2022 | Rank <= 5,
         Season != "Winter") %>% 
  arrange(Season) %>% 
  gt() %>%
    tab_header(
    title = md(paste0("**Warmest Seasons on Record**")), 
    subtitle = paste("Temperature Unit: Celsius")) %>%
  gt_highlight_rows(rows = Year == 2022, font_weight = "bold") %>% 
  tab_source_note(
    source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Data for the current season is controlled across years to only use the days reached in the current year.*"))
Warmest Seasons on Record
Temperature Unit: Celsius
Year Season Average Temperature Temperature Anomaly Seasonal Low Seasonal High Percent HW Rank
2012 Spring 7.97 2.15 5.46 12.90 100.00 1
2021 Spring 7.42 1.60 5.09 11.77 100.00 2
2022 Spring 7.41 1.59 5.65 11.92 100.00 3
2016 Spring 7.05 1.23 5.35 11.31 52.17 4
2013 Spring 7.04 1.22 4.68 11.33 71.74 5
2021 Summer 17.14 2.70 11.62 20.49 100.00 1
2012 Summer 16.80 2.39 11.14 20.24 95.29 2
2020 Summer 16.55 2.13 10.11 20.76 78.82 3
2022 Summer 16.53 2.09 11.83 20.43 94.19 4
2018 Summer 16.26 1.82 10.85 20.32 65.12 5
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Data for the current season is controlled across years to only use the days reached in the current year.

2022 Temperatures

Temperatures for the current year can be seen against the same thresholds with the following plot:

Code 
# Make Plot
hw_temp_p <- year_hw_temps(year_hw_dat = this_yr_hw, temp_units = "C") + theme_gmri() + theme(legend.position = "bottom") + labs(color = "")

# Show Figure
hw_temp_p

2022 Anomalies

Changing from absolute temperatures to anomalies reveals the degree to which this year is departure from what we would have expected from the climatology of 1982-2011’s temperatures.

Code 
# Show figure
hw_anom_p <- year_hw_anoms(year_hw_dat = this_yr_hw, temp_units = "C")+ theme_gmri() + theme(legend.position = "bottom") + labs(color = "")
hw_anom_p

Record Heatwave (Interactive)

Daily temperatures for the last 365 days are displayed below with reference to how they fall against the marine heatwave and cold spell thresholds.

Code 
# # Option 1: Grab last 365 days
# 
# # Grab data from the most recent year through present day to plot
# last_year <- Sys.Date() - 365 #1 year from current date
# last_yr_heatwaves <- region_hw %>% filter(time >= last_year)
# last_yr <- year(last_yr)


# # Option 2: Grab last Full Year:
# last_year <- Sys.Date() - 365 #1 year from current date
# 
# # wind it back to first day of the year
# last_year <- last_year - yday(last_year) + 1 
# 
# # Filter out:
# last_yr_heatwaves <- region_hw %>% filter(year(time) == year(last_year))
# last_yr <- year(last_yr)


# Option 3: Last complete year through now
last_year <- 2021
last_yr_heatwaves <-  region_hw %>% filter(year(time) >= last_year)


# Plot the interactive timeseries
last_yr_heatwaves %>% 
  plotly_mhw_plots()

3: Shifting Baselines

In 2021 NOAA is transitioning standard climatologies from the 30-year period of 1982-2011 to a new period spanning 1992-2020. Changes in climate regimes often does not result in a uniform upward or downward change that is consistent throughout the year.

The plot below shows just how both the average temperature, as well as the annual highs and lows have shifted. When looking specifically at Gulf of Maine here is how the expected temperature for each day of the year has shifted.

From this we can see that the Fall temperatures have risen more than those of the spring. There is also a large change in where the threshold for Marine Heatwave events sits, a consequence of exceptionally warm Fall temperatures becoming more common.

Code 
# Run heatwave detection using new climate period
heatwaves_91 <- pull_heatwave_events(region_timeseries, 
                                     threshold = 90, 
                                     clim_ref_period = c("1991-01-01", "2020-12-31"))  %>% 
  supplement_hw_data() %>% 
  filter(doy != 366)  %>% 
  mutate(year = year(time),
         yday = yday(time),
         yr_rev = factor(year),
         yr_rev = fct_rev(yr_rev),
         flat_date = as.Date("2000-01-01") + yday - 1)



# Subtract old from the new
heatwaves_91 <- heatwaves_91 %>% 
  mutate(clim_shift = seas - region_hw$seas,
         upper_shift = mhw_thresh - region_hw$mhw_thresh,
         lower_shift = mcs_thresh - region_hw$mcs_thresh)


# Make arrows where we want to point at things:
arrows <- 
  tibble(
    x1 = as.Date(c("2000-07-15")),
    x2 = as.Date(c("2000-08-28")),
    y1 = c(1.25), 
    y2 = c(1.15)
  )

# Plot the differences
heatwaves_91 %>% 
  filter(time >= last_year) %>% 
  mutate(year = year(time),
         yday = yday(time)) %>% 
  distinct(flat_date, .keep_all = T) %>%        
  ggplot(aes(x = flat_date)) +
  geom_line(aes(y = clim_shift, color = "Mean Temperature Shift")) + 
  geom_line(aes(y = upper_shift, color = "MHW Threshold Change")) + 
  geom_line(aes(y = lower_shift, color = "MCS Threshold Change")) + 
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(0.08, "inch")), size = 0.5,
    color = "gray20", curvature = -0.3) +
  annotate("text", x = as.Date("2000-06-01"), y = 1.225, label = "Fall Extremes\nMore Frequent") +
  labs(x = "", 
       y = "Shift in Expected Temperature \u00b0C",
       color = "") + 
  theme_gmri() +
  theme(legend.position = "bottom") +
  scale_color_gmri() +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))

 

A work by Adam A. Kemberling

Akemberling@gmri.org